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ABSTRACT 

We present the analysis of a suite of simulations of a Virgo mass galaxy clus- 
ter. Undertaken within the framework of standard cold dark matter cosmology, 
these simulations were performed at differing resolutions and with increasingly 
complex physical processes, with the goal of identifying the effects of each on 
the evolution of the cluster. We focus on the cluster at the present epoch and 
examine properties including the radial distributions of density, temperature, en- 
tropy and velocity. We also map 'observable' projected properties such as the 
surface mass density. X-ray surface brightness and Sunyaev-Zel'dovich signature. 
We identify significant differences between the simulations, which highlights the 
need for caution when comparing numerical simulations to observations of galaxy 
clusters. While resolution affects the inner density profile in dark matter simula- 
tions, the addition of a gaseous component, especially one that cools and forms 
stars, affects the entire cluster. For example, in simulations with gas dynamics 
but no cooling, improving the gravitational force resolution from 200 kpc to 14 
kpc increases the X-ray luminosity and emission-weighted temperature by fac- 
tors of 2.9 and 1.6, respectively, and it changes the form of the X-ray surface 
brightness and temperature profiles. At the higher resolution, a simulation that 
includes cooling and star formation converts 30% of the cluster baryons into 
stars and produces a massive central galaxy that substantially alters the cluster 
potential well. This cluster has 20% higher X-ray luminosity and 30% higher 
emission-weighted temperature than the corresponding cluster in the no-cooling 
simulation. Its properties are reasonably close to those of observed XD clusters, 
with conversion of cooled gas into stars greatly reducing the observational con- 
flicts found by Suginohara & Ostriker in simulations with cooling but no star 
formation. We conclude that both resolution and included physical processes 
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play an important role in simulating the formation and evolution of galaxy clus- 
ters. Therefore, physical inferences drawn from simulations that do not include 
a gaseous component that can cool and form stars present a poor representation 
of reality. 

Subject headings: Galaxy Clusters; Numerical Simulations; Star Formation 

1. Introduction 

With masses exceeding 10-^^ M©, clusters of galaxies represent the largest gravitationally 
bound objects in the Universe. For a number of decades, these systems have been the 
focus of intense observational and theoretical study in an effort to understand the 
evolving inter-relationships between the hot diffuse intracluster gas, the galaxies, and 
the inferred dark matter component that make up the clusters. In recent years, rapid 
advances in computer technology and available computational power, coupled with the 
development of increasingly sophisticated numerical simulation codes, have offered the 
possibility of studying the formation and evolution of cosmological structures like clusters 
in unprecedented detail. 

The developmental history of the increase in resolution and physical complexity of 
the numerical simulations has closely followed the growth of computational power. Early 
efforts considered simple systems consisting only of dark matter (e.g. Peebles 1970); the 
evolution of such systems depends solely on the action of gravity. Since the number of 
particles involved was relatively small, the gravitational forces were calculated directly. To 
enhance the degree of realism, novel techniques (such as particle-mesh and tree-codes) were 
developed so that simulations with an increasingly larger number of particles would be 
computationally feasible. 
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More recent simulations have also included a baryonic component. This was made 
possible using techniques like smoothed particle hydrodynamics (SPH), to study jointly 
the distributions of gas and dark matter within the cluster environment. Unlike the dark 
matter, the evolution of this baryonic material depends not only on gravity but also on 
processes such as cooling via radiative and Compton processes, heating via shocks and 
the deposition of energy by supernova explosions, the conversion of gas into stars, etc. 
Even in simulations where the baryonic component is treated as a simple adiabatic fluid, 
one still needs to consider the impact of processes like compressional and shock heating. 
Such heating modifies the subsequent motions of baryonic material within the cluster 
potential. The effects of these non-gravitational processes can be complex, and they can 
interact in unexpected ways with the effects of finite numerical resolution. The physics of 
star formation and supernova feedback is poorly understood; one can introduce plausible 
"recipes" for converting cooled gas into stars, but until these recipes are better constrained 
by observations they remain a source of uncertainty in numerical calculations. 

How much do the details of the resultant simulated clusters depend on the characteristics 
of the simulation? Specifically, how do properties such as the mass and force resolution and 
the degree to which non-gravitational physical processes are included affect the results? In 
this paper we address these questions by considering a family of simulations modeling the 
formation of a Virgo-sized cluster evolved with differing resolutions and including different 
physical processes. A complementary study was recently presented by Bryan and Norman 
(1998); while they considered a sample of clusters, their simulations were undertaken at 
lower resolution and did not include the cooling of gas and subsequent star formation. 

In §2 we briefiy review the results from recent simulations. In §3 we describe our 
methods. The three-dimensional radial profiles of a number of physical properties of 
the simulated clusters are compared in §4, more specifically the total and dark matter 
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distributions in §4.1, the circular velocity in §4.2, the gas density profile and baryon fraction 
in §4.3, the cooling timescales in §4.4, the gas temperature distribution in §4.5 and the 
entropy in §4.6. Projected, two-dimensional distributions, akin to observable properties 
of galaxy cluster, are presented in §5. These are the total and dark matter distributions 
in §5.1, X-ray surface brightness in §5.2, X-ray temperature distributions in §5.3 and the 
Sunyaev-Zel'dovich effect in §5.4. We discuss the implications of this study in §6, with §6.1 
focusing of the effect of resolution, §6.2 on the addition of a hydrodynamic component and 
§6.3 the effect of allowing the gas to cool and and form stars. The conclusions of our study 
are presented in §7. 



2. Past Simulations 

Many simulations of cluster formation have been undertaken over the last three 
decades. While the earliest had only a dark matter component, many have included a gas 
component and, in more recent work, have incorporated the physics of cooling and star 
formation within the gaseous environment. To exhaustively summarize the details of all the 
previous simulations is beyond the scope of this paper. For comparison purposes, however, 
we present the details of a representative sample of recent cluster formation simulations 
in Table 1. In addition to the physical details of the simulations, such as the final cluster 
mass and whether or not the gas was allowed to cool and form stars, we also present the 
resolution limits of the various components. This includes the spatial resolution, which we 
define to be twice the gravitational softening length, and the dark matter and gas mass 
resolution, which we define to be 32 times the mass of the corresponding particle. Wherever 
possible we convert the gravitational softening length into its equivalent Plummer softening 
value. Most simulations have a spatial resolution that is constant in comoving coordinates. 
We append a "p" to those that have a resolution constant in physical coordinates. For 
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the Eulerian simulations we quote the gas mass resolution at the mean density. The gas 
mass resolution for a virialized object is the number in parentheses. We also define a virial 
mass limit to be the mass of a virialized object, an object with a mean overdensity greater 
than 178, that has a radius that is twice the spatial resolution scale. This should be the 
minimum resolvable mass of a collapsed object. All values are at z — 0. As one can see 
from this table, the various resolutions used span several orders of magnitude. 

As noted in Suginohara & Ostriker (1998), these various methods reproduce the 
generic features of observed galaxy clusters, such as their core radii, suggesting that 
the overall picture of clusters forming from gravitationally driven merging and collapse 
of subcomponents provides a reasonable description for structure formation. While we 
consider this encouraging, we also feel that, owing to various numerical effects and the 
exclusion of various known physical components and processes, these simulations probably 
do not provide a complete description of the formation and evolution of galaxy clusters. 

It has been argued that cooling should be important only on scales smaller than the 
cooling radius; gas within the cooling radius can cool within a Hubble time, leaving gas 
outside the cooling radius unaffected. However, the loss of central pressure support owing 
to this cooling will influence the hydrodynamic evolution of the entire cluster and cooling 
during the early stages of cluster formation can remove gas from the eventual intracluster 
medium by converting it into stars. For galaxy clusters the cooling radius ranges from 100 
to 200 kpc depending on the size of the cluster. As we show below, cooling can greatly 
affect some physical properties of the cluster all the way out to the virial radius. 

Similarly, one might expect the limited spatial resolution to affect only scales smaller 
than the spatial resolutions listed in Table 1. Since most of the cluster simulations resolve 
the scales of observed cluster cores, one might perhaps conclude that all the simulations 
have sufficient resolution. Although this may be true for the total mass density distribution. 
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finite resolution effects modulate the potential out to several times the spatial resolution 
scale. Further, since clusters form hierarchically, one might expect that it is also necessary 
to resolve the subunits that eventually merge to form the cluster and that not resolving the 
subunits will influence even larger scales (the relevant comparable resolution is the fraction 
of the virial radius) . As we will show, some physical quantities are affected on larger scales 
than those listed in Table 1, sometimes approaching the virial radius. 

Judging by Table 1, the mass resolutions all seem adequate to resolve galaxy clusters. 
But once again, it may be necessary to resolve the substructures that merge to form the 
galaxy clusters. Most of the simulations appear to have mass resolutions sufficient to resolve 
several levels down the hierarchy. When cooling is included, however, it may be necessary 
to resolve galactic scales. 

3. Methods and Models 

For the study presented here, most of the simulations were performed with TreeSPH 
[Hernquist & Katz (1989), hereafter HK; Katz et al. (1996), hereafter KWH], a code 
that unites smoothed particle hydrodynamics [SPH; Lucy (1977); Gingold & Monaghan 
(1977)] with the hierarchical tree method for computing gravitational forces (Barnes & 
Hut 1986; Hernquist 1987). Dark matter, stars, and gas are all represented by particles; 
coUisionless material is influenced only by gravity, while gas is subject to gravitational forces, 
pressure gradients, and shocks. Because it uses a Lagrangian hydrodynamics algorithm 
and individual particle time steps, TreeSPH can perform simulations with the enormous 
dynamic range needed to study galaxy formation in a cosmological context. In SPH, gas 
properties are computed by averaging or "smoothing" over a fixed number of neighboring 
particles; 32 in the calculations here. When matter is distributed homogeneously, all 
particles have similar smoothing volumes. However, smoothing lengths in TreeSPH are 
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allowed to decrease in collapsing regions, in proportion to the interparticle separation, thus 
increasing the spatial resolution in precisely those regions where a high dynamic range is 
needed. In underdense regions, the smoothing lengths are larger, but this is physically 
reasonable because the gas distribution is smoother in these regions, requiring fewer 
particles for an accurate representation. To enable it to perform cosmological simulations 
TreeSPH includes periodic boundary conditions, comoving coordinates, radiative cooling, 
and an algorithm for star formation that turns cold, dense gas into coUisionless particles 
and returns supernova feedback energy to the surrounding medium. 

To maintain accuracy in the tree force calculation we use an opening angle criterion 
of of ^ = 0.7 (Hernquist & Katz 1989). We use a spline kernel for the gravitational force 
softening; this has an advantage over other forms for the softening (e.g. Plummer); forces 
become exactly Newtonian beyond twice the softening length. TreeSPH allows particles 
to have individual time steps according to their physical state, so that the pace of the 
overall computation is not driven by the small fraction of particles requiring the smallest 
time steps. To ensure accurate integrations we use the timestep criteria described in Quinn 
et al. (1998) and KWH for each particle with a Courant factor of 0.3 and eg^av — 0.4. 
To increase efficiency, we never allow the gas smoothing length to drop below 1/4 times 
the gravitational softening length. In the hydro dynamical calculations described here, the 
largest allowed timestep is 6.5 x 10^ years. The smallest timestep was 32 times smaller, or 
2.03 X 10^ years. 

To accurately simulate the formation of a galaxy cluster it is important to follow the 
evolution of the gravitational tidal field out to a large radius, while retaining high resolution 
in regions that eventually constitute the cluster. The global tidal field can affect the 
evolution of matter flowing into the cluster. To accomplish these goals using a minimum 
number of particles, hence making the problem computationally tractable, we use the 
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foUowing procedure. 

First, we used a Gelb & Bertschinger (1994) simulation of a large, uniform, 
periodic volume to create a catalog of galaxy clusters for further study. The simulation 
tracked the evolution of 144^ dark matter particles within a periodic cube 100 Mpc 
(Hq = 50kms~^ Mpc~^; h = 0.5 throughout this paper) on a side in a flat Einstein-de 
Sitter (Q = 1) universe. The mass of each dark matter particle in the simulation is 
Mdp = 2.3 X 10^° Mq. The dark matter particles were imprinted with perturbations 
described by the standard CDM power spectrum, and evolved to the present using a P^M 
algorithm (Gelb 1992) with a timestep of 1.9 x 10^ years and a (Plummer) gravitational 
softening length of 32.5h~^ comoving kpc. The amplitude of the initial perturbations is 
such that the square root of the variance in 8h~^ Mpc spheres (cg) is 0.7 when linearly 
extrapolated to 2; = 0. We selected one cluster, a relatively isolated system with a circular 
velocity of about 1000 km s~^, representative of a candidate Virgo-like cluster, for more 
detailed study employing a hierarchical mass grid for higher resolution (Katz & White 
1993). 

The Gelb & Bertschinger (1994) simulation is used to identify the particles that 
constitute the final cluster at the present day. The new simulation volume is centered on 
the eventual position of the cluster. The mass resolution in the finest hierarchical grid, a 
spherical region of radius 20 Mpc about the center of the simulation volume, is the same 
as that in the original simulation (i.e. Mjp = 2.3 x 10^° Mq), but it becomes successively 
coarser at larger distances from the cluster center. In this way, only 117,000 particles 
(100,000 in the central region and 17,000 outside) are needed to model the entire 100 Mpc 
cube, compared to over 3,000,000 if it were modeled at uniform resolution. We use the same 
realization of the power spectrum as before and evolve to the present using TreeSPH (no 
gas component) with a gravitational softening length of 35h~^ comoving kpc (equivalent 
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Plummer softening) and a maximum timestep of 5.2 xlO^ years. We use this simulation to 
initialize the final set of five simulations where the highest resolution regions are limited 
to only those containing particles that fall within the virial radius of the cluster at z = 
and their nearest neighbors. The details of these simulations are presented below and a 
summary of their resolution properties is presented in Table 2. 

The first of the five final simulations was a pure dark matter simulation, rerun with 
the same mass resolution in the high resolution regions as in the previous cases but with 
a gravitational softening length of 7h~^ comoving kpc (equivalent Plummer), a factor of 
five improvement in the spatial resolution. This simulation hereafter will be referred to as 
Dark. The total number of particles in this simulation is 68,000 with 44,000 high resolution 
particles and 24,000 particles of greater mass. 

We perform three further simulations at the same mass resolution as Dark but with 
the addition of a hydrodynamical gas component in the high resolution regions. Of these, 
two had the same gravitational force resolution and hence the same spatial resolution 
as Dark and the third had a gravitational softening length of 100h~^ comoving kpc 
(equivalent Plummer). Since we restrict the gas smoothing length to be greater than 1/4 
the gravitational softening length, this also has the effect of artificially reducing the spatial 
resolution of the gas. This low resolution SPH simulation and one of the two high resolution 
SPH simulations are evolved without supernova heating, radiative or Compton cooling, or 
star formation. We will refer to these two simulations as Lowres and Adiabatic, respectively. 
By comparing the results of these two simulations, we expect to assess the effects of spatial 
resolution. Note that we conform to the conventional misuse of terminology in this field 
by equating "non-radiative" gas dynamics with "adiabatic" gas dynamics, even though the 
gas in our adiabatic simulation is subject to shock heating and therefore does not evolve 
isentropically. 
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In the other high resolution hydrodynamic simulation (Cool+SF), the gas is allowed 
to undergo cooling via a number of mechanisms, including coUisional excitation and 
ionization, recombination, free-free emission and inverse Compton cooling off the microwave 
background. Coohng allows the gas to collapse, where possible, into cold, dense knots. In 
"real" systems, gas in such knots would form stars. Unhke gas, a stellar component is 
coUisionless and hence behaves very differently during collisions and mergers. Additionally, 
stars can inject energy into the surrounding medium via winds and supernova explosions. 
The fraction of the baryonic component locked up in stars can affect both the subsequent 
dynamical and thermal evolution of the system. Consequently, simulations that allow the 
gas to cool should also include a prescription for converting gas into stars; the Cool+SF 
simulation includes such a prescription for turning cold, dense gas into coUisionless "star" 
particles. The technique and its computational implementation are described in detail by 
KWH. In brief, gas becomes "eligible" to form stars if it has a physical density corresponding 
to nu > 0.1 cm~^ and an overdensity p/p > 56.7 (equivalent to that at the virial radius 
of an isothermal sphere). The gas must also reside in a convergent flow and be locally 
Jeans unstable, although the density criteria themselves are usually sufficient to ensure this. 
Eligible gas is converted to stars at a rate dlnpg/dt = —c^/tg, where tg is the maximum of 
the dynamical time and the cooling time. We use Cyr = 0.1 for the simulations here, but 
KWH show that the simulated galaxy population is insensitive to an order-of-magnitude 
change in c^, basically because the star formation rate is forced into approximate balance 
with the rate at which gas cools and condenses out of the hot halo. When star formation 
occurs, supernova heating is added to the surrounding gas assuming a standard IMF from 
0.1 to IOOMq and that stars above SMq become supernovae. Each supernova adds 10^^ 
ergs of thermal energy to the system. When a gas particle first experiences star formation, 
a fraction of its mass becomes stellar and the particle temporarily assumes a dual role as a 
"gas-star" particle, with its contribution to gas dynamical quantities dependent on its gas 
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mass alone. A gas particle is converted to a purely coUisionless star particle when its gas 
mass falls below 5% of its initial value. As discussed in KWH, the use of hybrid gas-star 
particles allows us to avoid the computational cost of introducing extra particles into the 
simulation but also avoids decreasing the local gas resolution in cooling regions, which 
would happen if we immediately converted gas particles into star particles. 

The initial gas particle distributions for the three SPH simulations are generated by 
taking the initial particle distribution for the Dark simulation and creating a gas particle for 
each high resolution dark particle. The gas particles are displaced from the dark particles 
by one-half of a grid spacing in each of the three directions. We assume a global baryonic 
fraction of Qb — 0.05 consistent with nucleosynthesis values (Walker et al. 1991). The 
resulting grid has 44,000 gas particles, 44,000 high resolution dark matter particles, and 
24,000 coarser resolution particles, for a total of 112,000 particles. 

Finally, we evolve a very high resolution dark matter only simulation (Hires) by 
splitting each particle in the Dark simulation that fell within the virial radius of the cluster 
at z — into 27 particles and extending the power spectrum down to smaller wavelengths. 
This simulation had a total of 1.3 million particles and was evolved using PKDGRAV 
(Stadel & Quinn 1998) with a gravitational softening length of 0.7h~^ comoving kpc 
(equivalent Plummer) and a dark matter particle mass of Mdp = 8.6 x 10^ Mq. 

All the results presented in this paper are a.t z = 0. The cluster virial radius, defined 
such that p{< Rvir) — 178pc, where pc is average density of a critical universe, is ~ 2Mpc 
in all the simulations. The cluster virial mass is 4.1 x 10^^ Mq and the circular velocity at 
the virial radius is ~ 1000 km s~^. 
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4. Cluster Properties: 3-D Distributions 

In this section we examined the physical properties of the simulated cluster averaged 
on spherical shells. During this section one must remember that the cluster is neither 
spherical nor even homogeneous on ellipsoidal shells. Interactively examining the structure 
of the cluster reveals sublumps, shocks and wakes whose properties are washed out in the 
spherical averages. For example, in the adiabatic simulation there is a sublump about 1 
Mpc from the cluster center that has a higher pressure and density but whose temperature 
is slightly cooler than the surrounding medium (lO^K vs. a few lO^K). Trailing behind this 
lump is a cooler "wake" (few lO^K) extending out to the virial radius, which is in rough 
pressure equilibrium with the surrounding hot gas. Ahead of the lump is a shock front with 
particles as hot as lO^K. Structures of this sort might conceivably be observable with the 
next generation of X-ray telescopes. Such substructure must also be kept in mind during 
the discussion when we refer to spherically averaged quantities. 

4.1. Total and Dark Matter Distributions 

In Figure 1, wc plot the mean radial profile of the total (dark + baryonic) mass density 
in the cluster, relative to the universal density, pc- The vertical bars along the x-axis 
denote the equivalent Plummer resolutions derived from the resolution lengths employed 
in the Adiabatic, Dark, Cool+SF and Lowres simulations. The heavier bars on top of the 
curves indicate the radii within which there are 32 particles; this can be considered the 
"mass resolution" limit of the simulation. At radii greater than ~ 300 kpc, the cluster mass 
density profiles in the different simulations are indistinguishable. At the virial radius they 
fall as p oc r~^. The effects of resolution and the included physical processes are apparent 
at radii less than ~ 300 kpc. The mass density profile in the Lowres simulation becomes 
increasingly shallower with decreasing radius, becoming nearly flat for r < 60 kpc. 
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The radial mass profile for the higher resolution simulation (Dark), does not appear 
to differ significantly from the Hires simulation except perhaps at r < 40kpc, suggesting 
that the total mass density profiles in these simulations are free from resolution effects for 
r > 40kpc. This agrees with the conclusions of Moore et al. (1998). The total mass density 
profiles in simulations that include the same physical processes appear to converge at radii 
greater than our defined spatial resolution limit. 

The inclusion of an adiabatic gas component does not significantly affect the total mass 
density profile on large scales. Differences do occur on small scales, however. The profile 
in the Adiabatic simulation slightly steepens within 60 kpc and then flattens substantially 
within 30 kpc, about the spatial resolution scale. The inclusion of cooling results in a 
significant change in the mass profile, with a rapid steepening within the central ~ 100 kpc, 
approximately 1.5 times the cooling radius in the Adiabatic simulation. 

The mass density profile can be approximated by 

P(^) ~ Tl \TT7-, , / / Nil. 6 (-'■) 
[r/rs) [l + {r/rs)\ 

where po ~ 5900pc and = 370 kpc. This profile has the same inner slope of o; = 1.4 found 
by Moore et al. (1998), although the slope in the outer region is slightly steeper with a ~ 3 
. [More recent high-resolution simulations of haloes of different masses all show an ~ r"^-^ 
inner profile, turning over to at large radii (Moore et al. 1998, 1999).] Equation 1 
matches the density profile of Navarro, Frenk & White (1997) at large radii but has a 
steeper central cusp. 

We plot the difference between the profile in equation (1) and the density distributions 
of our simulations in Figure 2. The parameterization of equation (1) accurately describes 
the Hires simulation, and it reasonably describes both the Dark and Adiabatic simulations 
except in the very central regions where these fiatten slightly. As indicated in Figure 1, 
both the Lowres and Cool+SF simulations have significant departures from the density 
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profile of equation (1). 

The total mass density profile for the cluster in the Cool+SF simulation, even though 
it has the same spatial resolution as the clusters in the Adiabatic and Dark simulations, 
continues to rise steeply towards the cluster center and appears to steepen even further 
within the central 20kpc. The density within the central 20kpc exceeds that in the 
Adiabatic and Dark clusters by a factor of 14. This behavior is caused by additional physical 
processes, in the form of cooling, star formation and feedback, included in this simulation 
and is caused by baryons collecting in the cluster center. This pooling of baryons in the 
cluster center deepens the potential well and modifies the dynamics in that region. 

In Figure 3, we show the radial density distribution for the dark matter component 
in each of the simulations. Since the Dark and Hires simulations are pure dark matter 
simulations, the curves appearing in this plot are identical to those in Figure 1. Interestingly, 
while all the simulations which include gas possess a very similar dark matter profile 
beyond ~ 30 kpc, within this radius it is apparent that the additional physical processes 
of gas cooling and star formation result in a radical modification in the form of the dark 
matter profile. This is not the case for those simulations where the gas was treated purely 
adiabatically, and the forms of the dark matter profiles for the Adiabatic and Lowres 
simulations are similar to those of their total mass distributions. 

In each simulation with gas, the baryonic component represents a minor fraction of the 
mass at all radii. The only exception is the very central regions in the Cool+SF simulation, 
where the baryonic fraction approaches Mf,/M 10, revealing the bayonic nature of the 
central regions of this simulation. While enhanced with respect to the other simulations, 
the dark matter in Cool+SF comprises only half the total mass density at 20 kpc and even 
less further in. 
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4.2. Circul£ir Velocity 

The circular velocity of a test particle at radius r in a spherical mass distribution is 
given by 

circ ^ \ J 

For an isothermal cluster, where the density falls as r~^, the circular velocity is constant 
with radius. Figure 4 illustrates the radial dependence of the circular velocity for the 
simulations presented in this paper. As in the previous plots, the vertical bars indicate the 
two spatial resolutions. The circular velocities derived from the various simulations are 
compared to that derived from the profile of Navarro, Frenk and White (1997) defined such 
that 

p[r) oc -- 



r (a + rY 



1 + log ( 1 + - 



(3) 



a-\-r J V a 

Such a profile is over-plotted in Figure 4, with Vmax ~ 1125 km s~^ at a radius of ~ 500 kpc. 
Unlike the total mass density profiles, which are identical at radii greater than ~ 300 kpc, 
the circular velocities for the different simulations only converge at r ^ 1 Mpc or slightly 
more than twice the nominal spatial resolution of the Lowres simulation. At the virial 
radius, the curves are very similar in form to Vnfw, decreasing as r~^ with increasing 
radius. 

For r < 1 Mpc, the circular velocity curves begin to differ from one another. This is 
very apparent for the Lowres simulation; Vc is maximum at ~ 1 Mpc and drops to ~ 85% 
of its peak value of ~1050km/s at the spatial resolution limit. The amplitudes of the other 
simulations all peak at ~1100km/s at 600kpc. Apart from slight variations in amplitude, 
the remaining curves track each other until ~ 200 kpc, where the circular velocity for the 
Cool+SF simulation stops declining as rapidly as the Adiabatic, Dark and Hires simulations. 
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The amplitudes and slopes of the latter three begin to diverge for r < 200 kpc. 

Between 100 kpc ~ r ~ 600 kpc the circular velocity of the CooZ-Z-^F simulation only 
decreases to ~ 95% of its peak value and then remains essentially constant until r ~ 20 kpc, 
where it begins to rise towards the cluster center While this behavior is expected given 
the nature of the total mass density profile, this region is within the resolution limit of 
the simulation. Such high circular velocities in the central 20 kpc of clusters are probably 
incompatible with the observed dynamics of brightest cluster galaxies. Conceivably better 
resolution would resolve this discrepancy, but it may instead point to a shortcoming in 
our input physics. The difference between this curve and that of the Adiabatic simulation 
illustrates the impact of baryon dissipation on the structure of the cluster potential. 

In sum, resolution affects the circular velocity out to typically 3 times the nominal 
resolution scale. The inclusion of an adiabatic gas component does not have much effect 
on circular velocities but a gas component that can cool and form stars causes the circular 
velocity to rise above pure dark matter calculations out to 200 kpc, about 3 times the 
cooling radius of the Adiabatic simulation. 

4.3. Gas Density Profile and Beiryon Fraction 

In Figure 5, we plot the radial gas density profile for the cluster in the Lowres, 
Adiabatic, and Cool+SF simulations. The two adiabatic simulations approximately converge 
at i? ~ 300 kpc, essentially the same as the convergence radius for the dark matter or total 
mass density profiles and corresponding to the nominal spatial resolution of the Lowres 
simulation. The gas density profile of the Cool+SF simulation, the one that includes cooling 
as well as star formation and its related feedback processes, docs not converge to the gas 
density profiles of the adiabatic simulations until the virial radius (~ 2 Mpc) . We have 
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already commented how these processes affect the structure of the very central regions of 
the cluster, with dissipation producing a strong density peak in the center. This density 
peak has a strong influence on the cluster's X-ray properties, as we discuss in § 5.2 and 
§6 below. Figure 5 shows that cooling and star formation affect the gas density profile 
throughout the cluster, partly by deepening the potential well, partly by converting a 
significant fraction of the gas into stars, and partly because these two effects in turn alter 
the hydrodynamic evolution of the remaining diffuse gas. 

Examining the profiles in more detail, one sees that all three gas density profiles have a 
slope of ~ 2.6 beyond the virial radius. At the virial radius the gas density profiles steepen, 
possibly indicating the presence of a weak shock. At 1 Mpc from the cluster center, both 
the Adiabatic and Lowres gas density profiles steepen to a slope of ~ 3.4. The steepening is 
less dramatic for the Cool+SF gas density profile which has a slope of ~ 2.7. Between 300 
kpc and the virial radius piowres > PAdiabaUc > PCooI+sf', to first order caused by the response 
of the gas to the different potentials. In the Lowres simulation the central potential is 
shallower and the gas distribution is more extended. At 300 kpc, piowres flattens to form a 
core and crosses the p Adiabatic profile, which continues to rise. At ~ 40 kpc, the gas density 
profiles of the Cool+SF and Adiabatic simulations also cross, owing to a combination of a 
rapid steepening of the Cool+SF density profile and a flattening of the Adiabatic proflle. 
In the Cool+SF case, gas cools to form a galactic component, much of it found in the 
cluster center in the form of stars. This effectively depletes the hot intracluster medium 
and forms a giant galaxy in the cluster center. The redistribution of the gas has important 
consequences for the baryon fraction in the cluster. The causes and implications of these 
processes will be discussed in more detail in §6. 

Figure 6 shows the fraction of the mass in baryons within a radius, R, normalized by 
r^b/il, the mean baryon density of the universe, which was 0.05 in the simulations, i.e. 
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we plot [Mb{< R)Qtot\/[Mtot{< R)^b\- We plot two curves for the CooZ-/- (SF simulation; 
the upper curve shows the distribution of the total baryonic component (in gaseous and 
stellar forms), and the lower curve shows the baryon fraction based only on the gaseous 
component. 

The baryon fraction curves from the Adiabatic and Cool+SF simulations behave 
differently. As expected, both curves asymptote to unity as the radius approaches and 
exceeds the virial radius. However, going from the virial radius to the cluster center, the 
baryon fraction in the Adiabatic simulation declines gradually until i? ~ 24 kpc and then 
levels off at ~ 20% of the universal value. In the Cool+SF simulation, the baryon (total) 
fraction curve rises steeply towards the cluster center. 

Looking at the gas fraction instead of the total baryon fraction in the Cool+SF 
simulation we find a very different result. Like the total baryon fraction in the Adiabatic 
simulation, the gas fraction in Cool+SF drops gently going from the virial radius towards 
the cluster center, leveling off at ~ 240 kpc. Moreover, the gas fraction in the cluster is 
significantly below unity at all radii and only climbs to 0.7 at twice the virial radius. The 
differences between the gas and the total baryon fraction curves in Cool+SF simulation 
indicates that 30% of the baryons in the cluster are locked up in stars and that the core of 
the cluster is, in fact, a stellar core. This is illustrated more graphically in Figure 7, where 
we plot the distribution of the individual baryonic components in the Cool+SF simulation; 
stars, hot gas (T > 10^-^K) and cold gas (T < 10^-^K). This plot is essentially unchanged 
when the boundary between hot and cold gas is reduced to lO'^'^K. Within 150kpc of the 
cluster center the baryons are dominated by stars, while beyond this radius, baryons exist 
mainly in the form of hot gas. Other than a small central peak and in galaxies , there is 
very little cold gas in the cluster; this has all been turned into stars. 



4.4. Cooling Timescales 



In Figure 8, we plot the cooling time as a function of distance from the cluster center. 
The cooling time is larger than the dynamical time at all radii in all three simulations. 
In the Lowres simulation none of the gas would cool within a Hubble time if cooling 
were suddenly included, something that might have led one to erroneously conclude 
that cooling is unimportant for the intracluster gas. Remember that in this and in the 
Adiahatic simulation cooling is not included during the evolution of the cluster. In the 
Adiabatic simulation, only within the central 70 kpc would the cooling time be shorter than 
the Hubble time (13 Gyrs) if cooling were suddenly included. Including cooling in this 
simulation would result in the eventual collapse of the gas within the central 70 kpc region. 

In the Cool+SF simulation, the only simulation that included cooling during the 
evolution of the cluster, the cooling radius is actually smaller, only 40 kpc, because the 
gas density has been reduced by the conversion of gas into stars. With our star formation 
prescription, the injection of energy from supernova explosions has no significant effect on 
the thermal properties of the intracluster gas; the cooling time scale in the central regions 
is very short and cold gas is readily turned into stars. Any energy that is injected into the 
gas rapidly escapes radiatively. Beyond the central region (> 900 kpc), sharp spikes are also 
apparent in Figure 8; these regions represents knots of cold, collapsed gas that correspond 
to galaxies within the cluster. 

4.5. Gas Temperature 

We plot the radial temperature profile of the gas averaged over a spherical shell in 
Figure 9. The temperature profiles for all the simulations are similar for radii greater than 
the virial radius, r ~ 2 Mpc, more than 5 times the nominal spatial resolution scale of the 
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Lowres simulation. At smaller radii they exhibit quite different behavior, however. The 
Lowres simulation has the shallowest rise towards the cluster center, leveling off at a value 
of ~2 keV in the central ~ 800 kpc. The cluster in this simulation has a large isothermal 
core. The central temperature hes slightly below the virial temperature of Tyir ~ 2.8 keV. 

The temperature profiles of the cluster in the Adiabatic and Cool+SF simulations 
begin to differ from one another within the central ~ 1 Mpc. The temperature profile for 
the Cool+SF simulation rises steeply towards the cluster center where the gas temperature 
exceeds 6 keV. The temperature profile for the Adiabatic simulation rises less steeply, 
reaching a maximum temperature of only 4 keV in the cluster center. 

The differences in the temperature profiles of the Lowres and Adiabatic simulations 
are caused by resolution effects, while the differences between the the Adiabatic and 
Cool+SF simulation result from the additional physical processes included in the Cool+SF 
simulation. Since one of these processes is cooling, we would have naively expected the gas 
in the Cool+SF simulation to be cooler towards the cluster center compared to the gas in 
the Adiabatic simulation. However, as we have already seen, the cooling time scale beyond 
40 kpc is longer than a Hubble time; consequently, coohng cannot directly affect the gas 
beyond the very central region. The indirect effect of cooling is to lower the gas density 
outside of the core and to raise the central concentration of baryons, and the deepening of 
the cluster potential due to the latter has the overall impact of raising the gas temperature 
throughout most of the cluster. Within 40 kpc, the gas temperature does plummet from 
the peak value of ~ 6.5 keV. We return to the impact of cooling on the cluster density and 
temperature profiles in §6. 

Comparing the gas temperature profiles with the gas density profiles, the gas density, 
in all three simulations, is related to the gas temperature by pgas T^ or equivalently, 
the pressure P oc pg{f for r > 800 kpc. In the central regions, the relationship approaches 
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Pgas oc T^/^ or P oc pgkf, the relationship for an isentropic gas. Between 800 kpc and 
40kpc the pressure-density relationship can be approximated by P oc p^ "^. These scaling 
relationships break down in the very central regions of the Cool+SF cluster. 

4.6. Gas and Dark Matter Entropy 

In the top panel of Figure 10, we plot radial profiles of a quantity proportional to the 
entropy per unit mass for an ideal gas, 

S = logio {^/P'£) , (4) 

where P = PgaskBTgas/A*niH is the gas pressure, Pgas is the gas density and //mn is the 
mean molecular mass of fully ionized gas {n — 0.6). In all three simulations, the mean 
entropy per unit mass of the gas at > 3 Mpc (i.e. well beyond the virial radius) is quite 
high — comparable to the entropy per unit mass inside the cluster. Why does the gas have 
high entropy at such radii? One possibility, suggested by the results of Cen et al. (1995) 
and Cen & Ostriker (1999), is that the high entropy of the gas outside the cluster may 
due to the gas having been heated to temperatures of 10^ — 10'' K by shocks associated 
with the formation of large-scale structures such as sheets and filaments. An alternative 
possibihty is that the high entropy of the gas outside the cluster is due to numerical effects, 
as we discuss further below. We will address the origin of gas entropy in more detail in a 
subsequent paper. Here, it suffices to note that the entropy per unit mass inside and outside 
the cluster differ by a very small amount. Therefore, accretion shocks must produce only a 
small amount of entropy, suggesting that the cluster accretion shock is weak. 

The gas entropy per unit mass in all three simulations differs within the virial radius. 
As we will discuss in §6, the total entropy of the gas in the three simulations is essentially 
the same. The profiles vary mostly because the gravitational potentials differ. That the 
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entropy profiles converge only at the virial radius highlights that differences in the potential, 
even if locahzed to the central 300 kpc of the cluster, can affect the gas distribution on 
scales of the entire cluster. The gas density profiles, as we already discussed, exhibit similar 
behavior. The entropy per unit mass decreases towards the cluster center, dropping from 
~ 7 at the virial radius to 5" ~ 5.5 at r = 40 kpc. In both the Adiabatic and the Lowres 
simulations, the intracluster gas forms a nearly isentropic core; in the Lowres simulation, 
the radial extent of the core is 300 kpc while radial extent of the core in the Adiabatic 
simulation is a factor of ~ 3 smaller. 

In the bottom panel of Figure 10, we plot radial profiles of a quantity characterizing 
the "coarse-grain" entropy per unit mass for the dark matter. 



where pdark is the local density of dark matter determined using the 64 nearest neighbours 
and (TiD is the one-dimensional velocity dispersion of those particles. 

On the whole, the dark matter entropy in all five simulations behaves similarly in 
the inner regions, declining from a value of ~ 5.5 at r = 0.8 Mpc to ~ 3.7 at r = 43 kpc. 
Over this same radial range, the dark matter radial density profile scales as Pdark oc r"^ 
for r > 0.1 Mpc and Pdark <^ for r < 0.1 Mpc. If we define Pdark = PdarkO"iD) then 
Pdark p}[^^.]^ for r > 0.1 Mpc, which is similar to the pressure- density relationship for the 
gas (Pgas oc Pg^). For r < 0.1 Mpc, however, the pressure-density relationship for the dark 
matter shows a fiattening trend (Pdarfc whereas that for the gas steepens to Pgas oc Pg^). The 
fiattening of the dark matter pressure-density relationship is due to the combination of two 
effects: One, the density profile is fiatter at r < 0.1 Mpc and two, the azimuthally-averaged 
dark matter velocity dispersion, which scales as r"'^-^ for r > 0.1 Mpc, turns over and 
decreases as r*^-^ towards the cluster center. 

Examining the dark matter entropy curves in detail, we note that the profile for the 
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Lowres simulation declines somewhat more gently towards the center than the corresponding 
profiles from the other simulations, a trend that appears to be directly correlated with the 
large force softening length in the Lowres simulation. The dark matter entropy profiles for 
the Dark, Adiabatic and Cool+SF also show a tendency to flatten on scales comparable to 
their respective force softening lengths. The entropy profile for the Hires simulation, which 
has much a smaller force softening length, shows no flattening. 

At distances greater than ~ 1 Mpc from the cluster center, the dark matter entropy 
profiles for four of the five simulations are similar. The only exception is the profile from 
the Hires simulation. Outside the cluster, the dark matter entropy in the Hires simulation 
is markedly lower than the entropy in the other simulations. We interpret the dark matter 
entropy outside the cluster as indicative of particle noise in the simulations. The particle 
noise is suppressed as the number of particles is increased, explaining the difference between 
results for the Hires and the other simulations. We note that since the number of dark 
matter particles and the number of gas particles in our simulations are equal, particle 
noise ought to contaminate both distributions to the same degree. In the case of the gas 
particles, however, artificial viscosity will convert the random motions due to noise into 
thermal energy when streams or shells cross, and it is possible that the high entropy of the 
gas beyond the cluster virial radius in the Adiabatic and Cool+SF simulations is due this 
type of particle noise heating. 

Also, the dark matter profiles at distances greater than ~ 1 Mpc have numerous 
troughs. These troughs are due to lower entropy dark matter in substructures. The Hires 
simulation exhibits the greatest degree of substucture and the corresponding entropy profile 
is highly serrated. In the three panels of Figure 11, we juxtapose the gas entropy profile and 
the dark matter entropy profiles. The offset between the two is mostly due to the difference 
in the dark matter and gas mass density. There is a high coincidence between the troughs 
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in the gas entropy and the dark matter entropy, again reinforcing the interpretation of the 
troughs in the dark matter distribution as being due to the presence of substructure. The 
number, the width and the depth of the troughs in the dark matter entropy profiles are 
directly related to the resolution of the simulations and the included physics. In comparing 
the Adiabatic and Lowres simulations, there are fewer troughs in the Lowres profile and the 
troughs are shallower and wider in keeping with the Lowres simulation having fewer, more 
"fluffy" clumps. 

Comparing the behaviour of the dark matter and gas entropy in the three simulations 
with gas, we find that, except for an offset, the two track each other very well on scales 
r > 0.2 Mpc. The offset between the two is due to the fact that the dark matter density is 
larger than that of the gas and the dark matter temperature is typically lower than the gas 
temperature by a factor of 0.8. On smaller scales, however, the dark matter entropy profile 
continues to decrease towards the cluster center whereas the gas entropy curve, at least in 
the Adiabatic and the Lowres simulations, flattens to form a nearly isentropic core. (The 
gas entropy profile in the Cool+SF simulation also shows signs of flattening in the inner 
0.1 Mpc but the steep decline in the gas entropy at the very center due to the cooling of 
the gas makes the comparison with the dark matter entropy less straightforward.) While 
this difference could be partly caused by caused by differences in numerical treatment of 
SPH and dark matter particles, we speculate that it is primarily a result of differences in 
the physical properties of the gas and dark matter components. In the case of the gas, 
short-range forces prevent colliding gas elements from interpenetrating, resulting in the 
formation of shocks. These shocks convert the kinetic energy of the gas into thermal energy; 
depending on the geometry of the shocks, the extent of acoustic waves generated by the 
shocks and the details of the artificial viscosity used in the simulations, the energy can 
be distributed over an extended region. Dark matter, on the other hand, is a coUisionless 
fluid and will interpenetrate. This interpenetration can result in dark matter elements with 
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very different velocities occupying tfie same small but finite spatial volume, resulting in 
a "mixed" system. The coarse-grain entropy, as we have defined it, is a measure of how 
well the system is mixed in this sense. The dark matter entropy profiles of our simulated 
clusters indicate that the dark matter in the central regions is not as well mixed as in the 
outer regions. This is also borne out by the decline in the dark matter velocity dispersion 
on scales smaller than 0.2 Mpc. The extent of mixing depends on the initial state of the 
system as well as on constraints imposed by integrals of motion (Tremaine, Henon and 
Lynden-Bell 1986) and therefore, is likely to be affected by the physical characteristics of 
the system, such as the angular momentum distribution of the dark particles, as well as by 
numerical effects, such as particle and force resolutions. 

5. Cluster Properties: Projected Distributions 

Simulations of galaxy clusters have a wealth of structure on a range of scales. 
Observationally, however, we are restricted to a two-dimensional view of these systems, in 
the form of X-ray luminosity and temperature distributions, and maps of the projected 
surface mass density derived from the analysis of gravitational lensing features. In the 
following section we present the radial distribution of projected quantities. For this, we 
consider a box with a length of 4 Mpc on a side, twice the virial radius, centered on the 
deepest point in the cluster potential. All quantities are first smoothed in three dimensions 
with an adaptive SPH-like kernel over the 32 nearest neighbors. This smoothed distribution 
is then projected along the required axis [see Katz & White (1993)]. 
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5.1. Total and Dark Matter Surface Mass Distributions 

The left-hand column of Figure 12 plots the radial profiles of the the projected total 
mass density and the right-hand column plots the radial surface mass density profiles for 
only the dark matter. Each panel, labeled X, Y and Z, plots projections along one of the 
simulation axes; the total mass within the box is 1.3 x 10^^ Mq. 

As expected from our analyses of the three-dimensional total and dark matter 
distributions in the different simulations, the projected total and dark matter mass 
distributions in the Lowres simulation appear similar. At large radii, the distributions fall 
off as R^^'^, which agrees with what one would have expected given the three-dimensional 
distributions. Towards the cluster center, the profiles flatten, eventually forming a core, 
which again is not surprising given the nature of the three-dimensional profiles. However, 
the size of the core depends on the viewing angle. In the X and Y projections, the core 
begins to form at ~ 300 kpc, consistent with what one would have expected given the 
three-dimensional profiles, but in the Z projection the core is smaller by a factor of ~ 3 and 
the central surface density exceeds that in the X projection by a factor of ~ 1.6 and that in 
the Y projection by a factor of ~ 2.6. It appears that this enhancement is not due to the 
fortuitous presence of a subclump of matter along the Z axis but rather it is caused by the 
central regions of the cluster being ellipsoidal rather than spherical in shape with the major 
axis aligned along the Z axis. 

The total and dark matter surface mass density profiles for the Adiabatic and Dark 
simulations are also very similar. There is a dependence of the core size and central peak 
projected mass density on the orientation of the cluster as seen in the Lowres simulation, 
although neither flatten enough to make a core in the central region. The projected mass 
density in the Adiabatic simulation has a slope of ~ 0.63 at radii less than 100 kpc in the X 
and Y projections while the projected mass density in the Dark simulation has a slightly 
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shallower slope of ~ 0.5. As in the Lowres simulation, the Z projection exhibits a markedly 
different behavior, with the Adiabatic simulation having a steeper slope of ~ 0.88 in the 
central regions and the Dark simulation having a slope of 0.66, both steeper than the slopes 
in the X and Y projections. The central projected mass density in the Z projection of 
the Adiabatic exceeds that of the other orientations by a factor of 1.7, while that of Dark 
simulation is 1.6 times higher. 

The projected surface mass density profile in the Hires simulation is virtually 
indistinguishable from its lower resolution counterpart, the Dark simulation. The projected 
dark matter and total mass distributions in the Cool+SF simulation are similar to the 
corresponding projections of the Adiabatic and Dark simulations for i? ^ 50 kpc. At smaller 
radii, however, the total surface mass density exhibits a sharp up-turn, with a slope of 
~ 1.15. This peak is not as prominent in the dark matter projections and is caused by the 
large central galaxy. Examining the profiles more closely, while both the total and dark 
matter surface density profiles in the X and Y projections are indistinguishable from the 
corresponding Adiabatic and Dark profiles beyond the central peak region, the projected 
mass density profile of the Cool+SF simulation in the Z projection lies below that of the 
Adiabatic and Dark profiles at clustercentric distances up to ~ 200 kpc. As in all the other 
simulations, the central surface density in the Z projection is greater than that in the other 
two projections, although the Z excess is only a factor of 1.25 in this case. 

Even though the mass within the virial radius is the same in each simulation, each 
exhibits different lensing signatures. It is possible for the Dark cluster to be super-critical 
to gravitational lensing, producing giant arcs, while the Lowres simulation lies below the 
critical mass threshold and produces no strong lensing features. The gravitational lensing 
characteristics of a galaxy cluster, such as the radial positions and widths of giant arcs 
and the distribution of weakly lensed images, depend on the core radius and central value 
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and slope of the projected mass distribution (e.g. Hammer 1991), not on the total cluster 
mass. A small, but very centrally condensed, galaxy cluster can produce a much stronger 
lensing signal than a more massive cluster that is very diffuse. Using coUisionless N-body 
simulations, Bartelmann et al. (1998) argue that the observed incidence of large arcs in 
clusters rules out most of the parameter space of CDM models, leaving only those models 
with low f^mattcr low Q\. They note that the effects of cD galaxies are a source of 
uncertainty in their conclusion, and our results suggest that the effects of radiative cooling 
and star formation could easily move some cosmological models from the "unacceptable" 
category to the "acceptable" category. 

Due to differences in the projected mass distributions the simulated galaxy clusters 
have a lensing signature that depends on the viewing angle. Therefore, the mass distribution 
inferred from an analysis of the gravitational lensing characteristics would also depend on 
the viewing angle. The X-ray profile, which we discuss in detail in the next subsection, does 
not depend strongly on the viewing angle and has about the same profile when viewed from 
all angles. This might partially explain the observed discrepancy between cluster masses 
derived from gravitational lensing techniques versus those derived from X-ray analyses 
(e.g. Allen 1998) and will be addressed more carefully in a forthcoming paper. 

5.2. X-Ray Surface Brightness 

As illustrated in Figure 9, the gas that permeates the galaxy cluster has a temperature 
of several keV, and such an ionized medium emits copious amounts of X-rays. Over 
the last two decades, observations of galaxy clusters with space-borne X-ray telescopes 
have provided important insights on the physical properties of intracluster gas. With the 
assumption of spherical symmetry and isothermality, a "fundamental plane" -like relation 
between the X-ray luminosity, temperature and mass (e.g. Edge & Stewart 1991) has been 
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used to probe the cluster potential, making X-ray analysis one of the primary tools for 
investigating clusters. 

The inclusion of a baryonic component into the simulated clusters allows us to similarly 
study such X-ray characteristics. The temperature and density of the gas is used, via 
a Raymond and Smith plasma code (Raymond & Smith 1977), to determine the X-ray 
emission at any location in the cluster. We project this emission along three orthogonal 
lines of sight, as we did with the mass distribution, to produce two-dimensional X-ray 
surface brightness maps using the SPH spline kernel smoothing [see Tsai et al. (1994) for 
details]. We assume the metal abundance of the gas is 0.3 Solar — a value considered typical 
for galaxy clusters (Mushotzky & Loewenstein 1997). 

We plot the X-ray surface brightness distributions for the simulated clusters in 
Figure 13. The left-hand panels plot the surface brightness distributions for a 2 — > 10 keV 
bandpass, similar to the ASCA Gas Imaging Spectrometers (GIS). We have not, however, 
filtered the X-ray surface brightness with any instrumental response functions and one 
should consider this to be a "raw X-ray" view of the cluster in this bandpass. This also 
allows a comparison to the X-ray luminosities in recently tabulated surveys of galaxy 
clusters [e.g. David et al. (1993)]. The right-hand panels plot the X-ray surface brightness 
distributions for a 0.1 2.4 keV bandpass similar to the ROSAT X-ray telescope filtered 
with the response function of the High Resolution Imager (HRI); these can be considered a 
true X-ray telescope view of the simulated cluster. 

For Cool+SF, Adiabatic and Lowres simulations, the X-ray luminosity of the clusters 
in the 2 lOkeV bandpass are 3.14, 2.58 and 0.88x10^^ erg s"-*^ respectively (these 
are also summarized in Table 3). In the same band pass the Virgo cluster has an 
observed X-ray luminosity of Lx = 1.16 x 10^^ erg s~^, with a bolometric luminosity of 
L^°' = 2.66 X 10^3 gj-g (David et al. 1993). The X-ray flux in this same bandpass is 
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within a factor of three of this observed value for all the simulated clusters. 

As with the total X-ray flux, the X-ray surface brightness profiles depend on the 
details of the simulation, but for each simulation these profiles have generic features that 
are independent of both the orientation and the energy band under consideration. Beyond 
1 Mpc, about half the virial radius, all the X-ray surface brightness profiles are similar. 
Within the inner ~ 300 kpc, the Lowres simulation profile rapidly turns over, forming a 
flat X-ray core. The X-ray surface brightness in the Adiabatic simulation continues to rise 
all the way into the center of the cluster, slowly turning over but never fully flattening to 
form an apparent core. The X-ray surface brightness profile of the Cool+SF simulation falls 
below the others for it! ^ 1 Mpc. Then, rising slowly, it crosses the X-ray surface brightness 
profile of the Lowres simulation inside its core. Within the central 50 kpc, about the cooling 
radius in the Adiabatic simulation, the profile of the Cool+SF simulation rapidly rises above 
the profile of the Adiabatic simulation. The cluster in the Cool+SF simulation has a bright 
X-ray core, a factor of 10 brighter than the central value in the Adiabatic simulation. 

The Cool+SF X-ray surface brightness profiles have a number of additional X-ray 
"spikes" caused by emission from cooling gas associated with galaxies. Such systems cannot 
form in the Adiabatic and Lowres simulations. However, superimposed upon their smoother 
profiles, one can identify additional bumps. These are associated with a prominent peak in 
the Cool+SF distribution. These deviations are caused by the interaction of the intracluster 
medium with a merging sub clump. 

In contrast to the projected mass density profiles, the viewing angle has very little 
effect on the shape or normalization of the X-ray profiles. The only exception is the location 
of the galaxy peaks in the Cool+SF profiles. The slopes of the X-ray profiles do depend on 
the energy band, however. 

To summarize, resolution affects X-ray surface brightness profiles within the nominal 
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spatial resolution scale, flattening them to form a core. The addition of cooling and star 
formation affects the X-ray surface brightness profiles out to one half the virial radius, both 
reducing their brightness and changing their shape to an approximate power law. Within 
the cooling radius, the profiles in the simulation with cooling and star formation rapidly 
rise, forming a bright X-ray core. Such sharp central peaks are observed in a substantial 
fraction of galaxy clusters and are interpreted as the observational consequence of a cluster 
"cooling flow" (see, e.g., Fabian 1994); whether or not this is the case in our simulated 
clusters will be addressed in §6. 

The X-ray surface brightness distributions of observed galaxy clusters are often fit by 
a simple parameterized model given by 



where r^ is the core radius and (3 controls the slope of the profile. Observed clusters have 
core radii of 100 — 600 kpc and /9 ~ 0.7 (Jones & Forman 1984). We fit this model to the 
2— >10 keV X-ray profiles in Figure 13. This model fits the Lowres simulation X-ray profiles 
with a core radius of 560 kpc and j3 — 1.4. 

Figure 14 shows fits of the /5-model to the three projected X-ray surface brightness 
profiles of the Cool+SF simulation (left-hand panel) and the Adiabatic simulation (right- 
hand panel). The presence of the central X-ray cusp (or bright galaxy) and the galactic 
spikes in the Cool+SF simulation complicates the fitting procedure. Following the analyses 
of X-ray emissions from XD clusters by Jones & Forman (1984), we excise the central cusp 
region. We also excise the galactic spikes. We use only those locations marked with a circled 
point in the fits. We fit two different /3-models to the simulations: one uses the full range 
of the X-ray profile and the second uses only a restricted range, from 50 kpc to 750 kpc, 
approximately the range used for observed clusters (Peres et al. 1998). We superimpose 
these fits over the profiles in Figure 14 with the dot-dashed hue corresponding to the fit over 
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the full range and the solid line to the fit over the restricted range. The /3-models provide a 
reasonable fit to the data, with the fit over the restricted range accurately reproducing the 
X-ray surface brightness within 0.5Mpc (neglecting the central cusp in the Cool+SF). The 
result of fitting to the restricted data range decreases the core radius and causes the slope 
to become shallower, i.e. /3 changes from ~ 1 to ~ 0.7. Observationally, the fitting of X-ray 
profiles find a /3 ~ 0.7, while measurements of f3 considering the gas temperature and cluster 
dynamics yield a value of ~ 1 [e.g. Jones & Forman (1999)], a situation often referred to as 
the /3-discrepancy. While several solutions to this problem have been put forward (Edge 
& Stewart 1991; Bahcall & Lubin 1994; Carlberg et al. 1996), the results found with our 
numerical simulation suggest that the situation may be resolved if X-ray profiles could be 
determined to larger radii [e.g. Navarro, Prenk & White (1995)]. The addition of coohng 
and star formation does not greatly affect the /3-model parameters when the central region 
is excluded. In fact, the model parameters are much more sensitive to the range over which 
the fit is performed than to the inclusion of cooling and star formation. 

5.3. X-Ray Temperature 

In addition to measuring the X-ray brightness, recent advances in X-ray detector 
technology have provided sateUite telescopes, such as ASCA, with increased spatial and 
spectral resolution, allowing the determination of X-ray emission weighted temperature 
maps of clusters (Markevitch et al. 1998). Such observations provide a more direct probe 
of the physical properties of the intracluster gas, revealing that the majority of clusters 
deviate substantially from isothermality and have strong temperature gradients, with many 
being a factor of two hotter in the central regions than at the virial radius. 

Since the X-ray emissivity as a function of position in the simulated cluster is known, 
it is straightforward to use the three-dimensional temperature information to determine 
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the projected X-ray temperature profiles, weighted with the X-ray emission in the chosen 
bandpass, for the simulated clusters. We plot these temperature profiles in Figure 15. 
The temperature profiles depend on the energy band because the emissivity in each band 
samples gas at different temperatures. 

In the 2 — > 10 keV band the Lowres simulation exhibits an almost flat, isothermal- like 
profile. The emission weighted temperature of ~ 2 keV is slightly lower than the virial 
temperature of Tvir ~ 2.8 keV. In the outer parts, r > 1 Mpc, the Cool+SF and Adiabatic 
temperature profiles gradually decline, falling by a factor of ~ 2 between 1 and 4Mpc, 
and hence deviate from an isothermal temperature profile. Within 1 Mpc, about half the 
virial radius, both temperature profiles begin to rise, peaking at ~ 2Tvir for the Cool+SF 
simulation and ~ 1.4Tvir for the Adiabatic simulation. The Cooi-Z-^F temperature profile 
falls to ~ l.STvir in the very central regions due to the presence of a small amount of cold 
gas. The profiles in the Adiabatic and Cool+SF simulations are similar to the profiles 
observed in galaxy clusters (Markevitch et al. 1998). The simulations have mean emission 
weighted temperatures of ~ 1.2Tvir for the Adiabatic simulation and ~ 1.6Tvir for Cool+SF 
simulation. However, in contrast to the spherically averaged 3-dimensional temperature 
profiles shown in Figure 9, projection causes the X-ray emission weighted temperature to 
have a less dramatic rise towards the center of the cluster. Also, the resolution of X-ray 
telescopes tends to be rather poor and much of the structure seen in Figure 15 would not 
be resolved, resulting in an apparently shallower gradient in the observed temperature 
distribution. 

In the ROSAT pass-band, the Lowres simulation no longer has an isothermal 
temperature profile but slowly rises to O.TTvir in the cluster center. Similarly, the Adiabatic 
simulation rises to 1.4Tvir, while the Cool+SF simulation rises to 2Tvir- In the very center, 
the temperature profile in the CooZ-Z-^F simulation drops sharply to 1.4Tvir, again due to the 
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presence of cold gas. This decline is more pronounced in the ROSAT band, since it is more 
sensitive to emission from colder gas (remember that the X-ray spectrum has been filtered 
with the instrument response of the ROSAT/HRI camera). At r ^ 1.4 Mpc, both the 
Cool+SF and Adiabatic simulations follow the temperature profile of the Lowres simulation, 
but within 1 Mpc both depart significantly, rising to peak temperature of ~ 1.4Tvii.. 

Finite numerical resolution affects the projected X-ray emission weighted gas 
temperature profile even on scales that are more than ~ 4 times the nominal spatial 
resolution scale. In the Lowres simulation this extends all the way out to the virial radius. 
Improving the spatial resolution by a factor of ~ 6 increases the mean emission weighted 
temperature of the cluster by almost 60%. Adding cooling and star formation increases 
the mean temperature by another 30%, to more than 1.5 times the virial temperature, 
and increases the central temperature by almost 50% over that found in the Adiabatic 
simulation. 



5.4. Sunyaev-Zel'dovich Decrement 

The bulk and thermal motions of the hot gas that pervades the cluster environment can 
affect the spectrum of the Cosmic Microwave Background (CMB) via Compton scattering. 
This induces an anisotropy in the observed CMB distribution, measured as a decrement in 
the Rayleigh-Jeans temperature of CMB radiation in the direction of the cluster (Sunyaev 
& Zeldovich 1972). Such Sunyaev-Zel'dovich decrements have been observed in several 
clusters at radio wavelengths [e.g. Myers et al. (1997)] , and more recently at sub-mm 
wavelengths (Lamarre et al. 1998), and are proving to be a useful probe of the intracluster 
medium, as well as offering clues to the values of cosmological parameters (Furuzawa et al. 
1998). 
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The inclusion of a baryonic component in the simulations allows us to investigate the 
Sunyaev-Zel'dovich decrement in the synthetic clusters. Figure 16 plots the combined local 
kinematics and thermal Sunyaev-Zel'dovich decrement, AT/T, as a function of radius (note 
that we do not consider the kinematic effect due to the bulk cluster motion). This quantity 
is independent of frequency in the Rayleigh- Jeans region of the spectrum. Each projection 
axis is the same as those in Figure 12. Only in the very outermost parts of the cluster, 
beyond the virial radius, does the decrement depend on the viewing angle, although even 
there the differences, and the Sunyaev-Zel'dovich signal, are small. In the central regions, 
R < 400 kpc, the individual profiles diverge, with the Lowres simulation turning over. Both 
the Cool+SF and Adiabatic simulations have decrements that continue to rise. Although 
the Cool+SF simulation appears shallower at larger radii, in the very center it peaks, with a 
central value similar to that of the more gently rising Adiabatic simulation. These features 
do not depend strongly on the cluster orientation. In the central regions of the cluster, 
the thermal component dominates the kinematic Sunyaev-Zel'dovich effect by an order of 
magntude. This is consistent with the small flow velocities (~ 50kms~^) seen in this region. 

6. Discussion 

It is apparent from the previous Sections that both numerical resolution and input 

physics can affect the evolution of a simulated cluster. Here, we consider the impact of 
these individually. 

6.1. Resolution 

Examining Figures 1 and 3 one can see that changing the force resolution affects the 
total mass and dark matter profiles only within about twice the spatial resolution scale. 
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The Lowres and Adiabatic simulations have identical profiles beyond that point. The mass 
profiles of the Lowres simulation show a significant flattening that begins on scales of twice 
the force resolution. The Adiabatic mass profiles continue to rise but also eventually flatten, 
again on a scale about twice the spatial resolution. The Dark simulation, consisting of 
only dark matter, has the same spatial resolution as the Adiabatic simulation. Compared 
to the Hires simulation, which is also a pure dark matter simulation, but with a spatial 
resolution ten times smaller than the Dark simulation, once again both the total mass 
and dark matter profiles of each simulation are the same until approximately twice the 
Dark simulation spatial resolution scale, where the Dark simulation mass profiles flatten 
slightly and depart from the Hires simulation mass profiles. Moore et al. (1998) also find a 
continual steepening of the mass profile as the spatial resolution is increased in their dark 
matter only simulations, but the mass profiles eventually converge as the spatial resolution 
scale became small enough (at fixed mass resolution). However, they also find that the 
mass profile becomes steeper as the mass resolution is increased, which accounts for some 
of the differences in the mass profiles of the Dark and Hires simulations. 

The degree to which particles can bind together gravitationally is limited by the 
hardness of gravitational interactions, resulting in deeper potentials for smaller softenings. 
But mass resolution also plays a part: the smaller the particle mass, the denser the central 
core can become, steepening the potential. As mentioned in Section 5.1, the steepness of 
the central core can affect the gravitational lensing characteristics of a cluster of galaxies 
[e.g. Bartelmann (1996)] and one might draw erroneous conclusions about the expected 
lensing signature from a model cluster if the mass profile has not converged. 

While resolution affects the form of the mass density profile on scales of twice the 
spatial resolution, the cluster potential, which depends on the integrated mass profile, is 
modified on much larger scales. Therefore properties that probe the potential, such as the 
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circular velocity and the temperature, are affected more globally. In Figure 4, the difference 
between the central mass profiles of the various simulations leads to different circular 
velocity curves on scales approaching half the virial radius. The distribution of gas within 
the cluster environment also depends on the form of the cluster potential. In Figure 9, the 
shallower cluster potential in the Lowres simulation results in less compressional heating 
than the higher resolution simulation whose potential is deeper, causing the temperature 
to be affected out to the virial radius. Finally, the total X-ray luminosity of the cluster 
also depends on the resolution. The higher resolution Adiabatic cluster has, as expected, a 
higher X-ray luminosity than the Lowres cluster. 

6.2. The Addition of Baryons 

Before discussing the detailed evolution of the baryonic material in these simulations, 
it is prudent to investigate whether or not the gas is artificially heated by numerical effects, 
specifically two-body interactions similar to those responsible for heating and subsequent 
core collapse in stellar systems. Such effects were recently investigated by Steinmetz & 
White (1997), whose analysis revealed that such heating effects can significantly modify the 
properties and hence evolution of the gaseous component. The magnitude of this effect 
depends mainly on the mass of the dark matter particles, although the relative gas particle 
to dark matter particle mass can mildly influence the cluster evolution. Steinmetz & White 
(1997) demonstrate that if the mass of the individual dark matter particles is below that 
of a critical mass, given by Merit ~ 5 x IO^VT^M©, where Tg is the temperature in units of 
lO^K, then spurious heating affects should be unimportant. In the high resolution regions 
of our simulations, where the temperature is several keV, Merit ~ 4 x lO^^M©; the masses 
employed in the simulation are Mgp = 1.2 x 10^ M0 and M^p = 2.2 x 10^^ Mq, respectively. 
We conclude, therefore, that our simulations should be free from spurious two-body gas-dark 
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matter heating effects. 

In comparison to the Dark simulation, including a baryonic component in the form of 
an adiabatic gas (Adiabatic) does not greatly affect the overall properties of the cluster. 
Looking at Figures 1 and 3, the total mass and the dark matter profiles both have a very 
similar form outside of a radius a few times the spatial resolution. Within this distance, the 
mass profiles of the Adiabatic simulation rise slightly more steeply than the Dark simulation, 
before flattening. The Dark profile continues to rise into the center of the cluster. As with 
the effects of resolution, this redistribution of matter in the cluster, and hence modification 
of the cluster potential, influences the evolution of the cluster on scales several times the 
spatial resolution. 

6.3. The Addition of Cooling and Star Formation 

The inclusion of additional physical processes in the form of cooling and star formation 
can radically affect the evolution and subsequent properties of a cluster. For example, 
consider the proflles of the total mass and dark matter distributions (Figures 1 and 3). 
The Adiabatic simulation flattens as it approaches the center, but the Cool+SF simulation 
continues to rise rapidly and has a central peak. Figure 7 indicates that this baryonic core 
is mainly stellar. 

How did the baryons pool into the center? Two possibilities present themselves, both 
of which may be occurring in the simulated cluster. First, the dense gas at the center of the 
cluster can dissipate its energy by radiative cooling, collapse to very high density, and form 
stars - in essence, a "coohng flow" creates a massive central galaxy. Second, dynamical 
friction and stripping of subclumps as they plunge into the cluster center can deposit stellar 
material into the cluster core (Babul & Katz 1993). Before the cluster itself assembles. 
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cooling allows gas to condense into the centers of "subunits," individual galaxies or galaxy 
groups. If the gas is cold and dense enough, it turns into stars. As a subunit merges into 
the cluster its less dense outer parts, comprised mostly of dark matter, are tidally stripped 
first and deposited further out in the cluster. Dynamical friction drags the dense core of 
the subunit closer to the cluster center, where eventually the tidal field is strong enough 
to disrupt it, depositing this baryon rich material deep in the cluster potential well. This 
process can cause segregation of baryons towards the center of the cluster even if the cooling 
times in the cluster are long. 

To fully differentiate between these possibilities one needs to trace the dynamical 
evolution of the baryonic material during the formation and subsequent evolution of the 
cluster. This is beyond the scope of this paper which focuses upon the end-point of the 
cluster simulations, but will be the subject of a forthcoming article in which the cluster 
properties will be traced through the earlier simulation outputs. We can, however, address 
the question of whether or not the Cool+SF simulation has a cooling fiow at the present 
time. We plot the gas mass accretion rate through spherical shells into the cluster, defined 
to be M = 47r/9pa5 (r)r^V^(r), normalized by the rate required for the cluster to accumulate 
its mass of gas within the virial radius (e.g. QeMyir) steadily over a Hubble time (i.e. 
M Ri 1 X lO^Mo/yr) in Figure 17. Each of the simulations has similar rates of accretion, 
with both the Cool+SF and Adiabatic simulations having essentially identical profiles into 
40 kpc. Since there is no cooling in the Adiabatic simulation, this implies that presently there 
is no significant cooling fiow in the Cool+SF simulation on a scale > 40 kpc. Considering 
Figure 13, it is within this radius that the X-ray profiles of the Coo/-/- 5"^ simulations possess 
an X-ray cusp; when coupled with the increase in the mass accretion in this region, the 
evidence points towards the Cool+SF simulation possessing a cooling fiow within 40 kpc. 
This scale is smaller than typical scale of the cooling region [150 kpc of the cooling flow 
clusters Peres et al. (1998)]. This cooling flow, however, has no signiflcant impact on the 
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inflow of mass on larger scales, suggesting that the latter is the result of recently accreted 
material being brought to a halt and redistributed with the cluster. 

We also examine the radial infall velocity of the gas in the three simulations. In 
Figure 18 we plot the ratio of the spherically averaged radial infall velocity to the dynamical 
velocity, defined as r/tdyn where tdyn = (Gptot)"^^^- The mean radial infall velocity is 
approximately 350kms~^ at distance of 4Mpc from the cluster center and scales roughly 
as \Vr\ oc r from 4Mpc to ~ 0.5 Mpc in all three simulations. Therefore, the characteristic 
time scale for the radially infalling flow over this range in radii is nearly constant at 
tflow ~ 10 Gyr. For r > 3 Mpc, the flow time is equal to the dynamical time, and the flow can 
be characterized as gravitational infall. Inward of 3 Mpc, pressure forces begin to affect the 
flow, and the flow velocity drops below the dynamical velocity. At r ~ 0.5 Mpc, the radial 
velocity of the flow is only ~ 50kms~^, and inward of 0.5 Mpc \Vr\ continues to decrease. 
In Figure 19 we plot the spherically averaged radial velocity of the gas divided by the 
local adiabatic sound speed, i.e. the Mach number. The flow undergoes a relatively rapid 
transition from subsonic to mildly supersonic flow in the region 2.5 Mpc < r < 3.5 Mpc, 
suggesting that a weak accretion shock is in this region. This places the shock at a slightly 
larger radius than that implied by the rise in the gas density profiles. The main differences 
between the results for the Lowres simulation and the other two simulations are caused by 
the smaller sound speed in the Lowres simulation. In the very central regions, the Lowres 
simulation appears to have a radial outflow. This feature is well below the spatial resolution 
of this simulation and is a numerical artifact. 

At ~ 1 Mpc the radial velocity of the dark matter component in each of the simulations 
has a strong dip, indicating a more rapid infall of material at this radius caused by a 
significant clump of material currently merging with the cluster. Closely examining the 
gas radial infall plots, the gas in both the Lowres and Adiabatic simulations has similarly 
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enhanced infall velocities at this radius, implying that the infalling clump contains both 
dark matter and gas. This feature, however, is not present in the gas of the Cool+SF 
simulation, as most of the gas in the merging clump has collapsed and turned into stars. 
Such infalling clumps, coupled with the action of dynamical friction, could represent a 
primary mode of transport of baryons into the central regions of the cluster. 

What causes the differences in the gas temperature profiles between the Cool+SF and 
Adiabatic simulations? We conjecture that two processes, either singly or jointly, could 
be responsible for these differences. First, supernova feedback could have heated the gas. 
The radial profile of the cumulative baryon fraction (Figure 6) shows that a non-negligible 
fraction of the baryons is in the form of stars. As a result of long local cooling time scales, 
the energy injected by the associated supernova explosions might not be radiated away 
and might heat the gas. The second possibility is an indirect impact of cooling. Cooling 
allows some of the gas to collapse and form stars, removing its contribution from the 
pressure support of the intracluster gas. Much of this collapsed, stellar component forms 
the dense, dynamically significant object in the center of the cluster. The gravity from this 
object draws in both the dark matter and the gas in the cluster; the drawing in of the 
former manifests itself as a steepening in the dark matter density profile. The gas would 
be similarly drawn in but unlike the dark matter would suffer compressional heating. Due 
to the long cooling time scale, the resultant thermal energy would not be radiated away. 
There is cooler gas in the central object, as well as in other objects that one would identify 
as galaxies, but the scale of these objects is much smaller than the scale of the cluster, and 
cold, dense gas turns into stars. 

The radial profiles of the entropy, S{r), for the Adiabatic and Coo/-/-<S'F simulations (see 
Figure 10) suggests that the major differences in the gas density and temperature profiles 
for the two clusters (outside the very central core region) are caused by compressional 
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heating of the gas being drawn adiabatically into the deep central potential rather than 
by supernova heating. It is as if a given mass (gas) shell at some radius Ri and entropy 
SAMg in the Adiabatic simulation has simply been moved to a radius R2 < Ri in the 
Cool+SF simulation. To verify that this is indeed the case, we determined the mean 
entropy per unit mass [/ S{r) dMg(r)/Mg(r < r-vir)] between the clustercentric radius R 
and the virial radius, as a function of the gas mass in the same volume. This quantity 
is unaffected by star formation. Under the above hypothesis, the results for the Cool+SF 
and Adiabatic simulations should be indistinguishable. As shown in Figure 20 this is 
indeed the case, and the results for the Lowres simulation are the same. Since the mean 
entropy per unit mass of the intracluster gas in all three simulations is essentially the same, 
supernova feedback, present only in Cool+SF simulation, cannot be significantly affecting 
the intracluster medium. The differences in the gas density, gas temperature, and specific 
entropy profiles are caused by the differences in the nature of the potential wells and by the 
gas responding adiabatically to establish hydrostatic equilibrium in these potentials. Hence 
one can understand the distribution of gas in the various simulations in terms of individual 
shells of gas behaving adiabatically. This does not imply, however, that the intracluster gas 
has a uniform entropy distribution (e.g. Figure 10). The entropy profile depends on the 
accretion history of the gas, including the strength of shocks. 

The gas distributions within the simulations presented here can be described by a 
polytropic equation of state with 7 1.1-1.2. Such a conclusion was also reached for 
a sample of nearby X-ray emitting clusters using ASCA observations (Markevitch et al. 
1998). In their study of cluster mass estimates using gravitational lensing and X-ray 
methods, Miralda-Escude & Babul (1995) found that the assumption that the intracluster 
gas is isothermal (as was thought at that time) results in a discrepancy between the two 
mass estimates, while the assumption that the gas was polytropic with 7 fa 1.1 resolves 
the discrepancy. The similarity between these observational results and our simulations is 
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encouraging. The poly tropic nature of the gas appears to be a robust result, since it occurs 
in all three simulations even though they differ in their resolution and in the included 
physical processes. We will further investigate the origin of such a distribution in our next 
paper, where we study the physics of cluster formation by examining the history of the dark 
matter and gas properties through the course of the simulation. 

Unlike simple gravitational physics and the hydrodynamics of an adiabatic gas, the 
behavior of a gas that cools and forms stars could depend on the prescription one uses 
to model these processes in a simulation. Although the prescription used in the Cool+SF 
simulation may not be a completely accurate description of star formation, the analysis 
presented here shows that cooling and star formation cannot be neglected. Consider, for 
example, the X-ray flux from a simulated cluster. The X-ray luminosity for the Cool+SF 
simulation is Lx ~ 3 x 10"^^ erg s~^, about three times that measured for the Virgo cluster of 
galaxies (which the simulation was chosen to represent) . This contrasts with the simulations 
of a 7 X 10^^ Mq cluster by Suginohara & Ostriker (1998), who find an X-ray flux of 
7.4 X lO^^erg s^^ with a temperature of only 3 keV. Clusters with this temperature are 
observed to have a much smaller flux, by more than a factor of 10 (David et al. 1993). 
While the Suginohara & Ostriker (1998) simulations and the CooZ-/-<S'F simulation presented 
here both allow gas to cool, theirs did not include star formation. Star formation turns 
cooler, dense gas into stars and removes it as a significant source of X-rays, and it appears 
that this additional physics goes a long way to repairing the discrepancy between predicted 
and observed X-ray properties found by Suginohara & Ostriker (1998). Therefore we do 
not require additional heating via conduction or an excessively high supernova rate to 
reproduce observed X-ray luminosities. However, while the reasonable X-ray properties of 
our simulated clusters are encouraging, it is also clear that these properties depend in detail 
on the resolution and physics of the simulation. In the Cool+SF simulation, the reasonable 
X-ray luminosity has been achieved through the removal of some of the X-ray emitting 



-46- 



gas by conversion into stars. The fraction of baryons in stellar form, however, is 30% as 
opposed to the observed ~10% and even then, Cool+SF cluster has the largest total X-ray 
luminosity. 

In previous simulations of clusters that have included cooling, but not the subsequent 
formation of stars, the pooling of gas into the cluster core has resulted in the formation of a 
bright central X-ray peak approximately 200 kpc in extent (Katz & White 1993; Suginohara 
& Ostriker 1998). As demonstrated here, the addition of star formation and supernova 
feedback into the baryonic component reduces the brightness and the size of the X-ray cusp 
but does not eliminate its formation. In our simulation, the X-ray cusp is the result of 
stellar baryons forming a central object massive enough to influence the dark matter profile 
significantly. On the basis of the X-ray surface brightness profile, the Cool+SF cluster most 
closely resembles the XD clusters studied extensively by Jones and Forman (1984). X-ray 
cusps are seen in ~ 70% of observed clusters, with the enhanced central emission interpreted 
as being due to a coohng flow (Peres et al. 1998). While the Coo/-/-<S'F simulation does 
exhibit a central X-ray cusp as well as evidence for a cooling flow in the central 40 kpc, the 
flow has no significant impact on the kinematics beyond 40 kpc. 

Even with the inclusion of cooling and star formation, simulated clusters of galaxies 
may not truly resemble 'real' galaxy clusters. However, this does not imply that the complex 
processes associated with baryonic matter can be neglected in numerical studies of galaxy 
clusters. If simulations with cooling and star formation do not reproduce the observed 
properties of galaxy clusters, it may be because a much more sophisticated treatment 
of cooling is required, both in physical modeling or numerical treatment, or it may be 
that the simulations neglect other important processes like magnetic fields permeating 
the intraclustcr medium or feedback from AGN. A similar conclusion was also reached by 
Suginohara & Ostriker (1998). 
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7. Conclusions 

• Moore et al. (1998) demonstrated that changing the resolution of a simulation results 
in a change in the slope of the density profile within the spatial resolution scale. We 
confirm this for the simulations presented here and similarly find that the central 
density profile is approximately p oc r~^^^. With the inclusion of cooling and star 
formation, however, the form of the mass profile is radically different, having a very 
steep central cusp. The combination of this deeper central potential well and the 
reduction in gas density in the cluster's outer parts (because of conversion to stars) 
affects many of the cluster's physical and observable properties. Also, we find that 
our lower resolution simulation does not properly capture the central density profile, 
illustrating that sufficient resolution is required to obtain convergent results. 

• The scaling of the three-dimensional gas temperature with gas density is similar to 
that deduced by Markevitch et al. (1998) for a sample of nearby X-ray clusters. This 
result appears to be independent of the details of the simulation. The gradient in 
the observed temperature profile depends on projection effects (as well as on the 
resolution of the X-ray telescope). 

• The X-ray surface brightness profile and total luminosity depend sensitively on both 
resolution and physics. In the simulation where gas cools and forms stars, a bright 
central X-ray cusp develops in the X-ray surface brightness profile. The total X-ray 
luminosity is reasonable for a Virgo-like cluster; however, the fraction of baryons in 
the form of stars is 30%, which is higher than typical fraction in observed clusters by 
a factor of ~ 3. Central X-ray cusps are observed in ~ 70% of galaxy clusters and are 
interpreted as being due to cooling flows (Peres et al. 1998). While a flow is apparent 
in the Cool+SF simulation, its scale (~ 40 kpc) is smaller than those typically seen in 
observed galaxy clusters (~ 150 kpc). Based on the X-ray surface brightness profile. 
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the Cool+SF cluster most closely resembles an XD cluster. 

• The fitting of /3-models, a common practice in cluster mass determinations from X-ray 
data, is sensitive to the radial range over which the fit is made. In our simulations, 
the best fitting model out to beyond the virial radius has ~ 1, but a fit over 
the restricted radial range typically observed yields [3 ~ 0.7. This suggests that 
the observed discrepancy between (5 measured from cluster X-ray surface brightness 
distributions and that determined from X-ray temperatures and cluster dynamics 
may be resolved if the fit to the observed X-ray surface brightness profile could be 
extended to large radii. 

• The gravitational lensing characteristics of a galaxy cluster, such as the radial positions 
and widths of giant arcs and the distribution of weakly lensed images, depend on the 
core radius, central value, and slope of the projected mass distribution (e.g. Hammer 
1991), rather than on the total cluster mass. We find, due to the non-spherical 
nature of the clusters in our simulations, that these properties also depend on cluster 
orientation with respect to an observer. Hence, the mass distribution inferred from an 
analysis of the gravitational lensing characteristics would also depend on the viewing 
angle. The X-ray profile, however, does not depend strongly on the viewing angle. 
This might help to partially explain the observed discrepancy between cluster masses 
derived from gravitational lensing techniques and those derived from X-ray analyses 
(e.g Allen 1998). The increased central concentration in simulations with radiative 
cooling and star formation may rescue some cosmological models that otherwise fail 
to reproduce the observed incidence of cluster arcs (cf. Bartelmann et al. 1998). 

• Due to computational limitations, simulations of cluster formation that have included 
a baryonic component have usually treated it as a simple adiabatic gas. Given that 
gas in clusters must have undergone cooling, collapse and processing into stars, is 
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this assumption of adiabaticity valid? Within our prescription both coohng and 
star formation proceed rapidly, resulting in very little cold gas pervading the cluster 
environment. The remaining hot gas continues to behave as an adiabatic fluid 
(Figure 20). However, ignoring these processes does not result in the same cluster 
evolution. Significant differences exist in all cluster properties for the Adiabatic and 
Cool+SF simulations. In the latter, as gas enters the forming cluster its density 
rises, so it cools efficiently and turns into stars. The reduction in pressure support 
in the central regions causes more gas to pool into the center and undergo cooling. 
This modification of pressure support, coupled with the increasing central density 
caused by the accumulation of stars, dramatically infiuences many of the cluster 
properties on scales out to the virial radius. It is important to note, therefore, that 
while the prescription used in the simulations presented here may not completely 
describe the cooling and subsequent star formation in galaxy clusters, the effects of 
these processes cannot be neglected if one wishes to draw physical inferences from 
numerical simulations. 
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Fig. 1. — The radial dependence of the density of all matter, both the dark and baryonic 
components, in each of the simulations presented in this paper. The bars along the x-axis 
indicate the spatial resolution for the medium and low resolution simulations. The equivalent 
bar for Hires lies off the left of the figure. The heavy bars that cross the individual curves 
indicate the radius within which each cluster has only 32 particles. 

Fig. 2. — The difference between the distributions in Figure 1 and the model given by 
Equation 1. 

Fig. 3. — As in Figure 1, but plotting the radial dependence of only the dark matter density 
in each of the simulations. Again, the bars along the x-axis indicate the softening lengths for 
the medium and low resolution simulations, while the heavy bars indicate the radius within 
which the cluster has only 32 particles. 

Fig. 4. — The radial dependence for the circular velocity in each of the cluster simulations. 
The bars along the x-axis indicate the spatial resolution for the medium and low resolution 
simulations. The thick solid line in this figure shows the circular velocity of a halo with an 
NFW profile. 

Fig. 5. — The radial dependence of the gas density in each of the simulations. The bars along 
the X-axis indicate the softening lengths for the medium and low resolution simulations and 
the heavy bars crossing the curves indicate the radius within which the simulation has only 
32 particles. 
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Fig. 6. — The radial dependence of the cumulative baryon fraction in each of the simulations. 
The two solid lines represent the Cool+SF simulation: the total baryonic content of stars 
and gas (top line) and only gas (bottom line). The bars along the upper x-axis indicate the 
spatial resolutions of the low and medium resolution simulations. The central cusp in the 
Lowres simulation is well within the 32 particle resolution limit. 

Fig. 7. — The distribution of the various baryonic components in the Cool+SF simulation. 
Here, the density is in units of the critical density and the black line represents the total 
baryonic mass. The central peak is comprised mainly of stars (dark grey line), although 
hot gas (lighter grey line) becomes the dominant baryonic component beyond ~ 130 kpc. 
Other than a peak at the very center of the cluster, there is very httle cold (T < 10^-^K) gas, 
represented by the hghtest grey hne, in the cluster. This has been converted into stars. 

Fig. 8. — The cooling time scale for the intracluster gas in the three simulations with baryonic 
component. Note the change of scale from logarithmic to linear at IMpc. The sharp spikes 
apparent in the Cool+SF simulation are knots of cold, collapsed gas. This gas is being 
processed into stars and represents galaxies. 

Fig. 9. — As in Figure 1 but plotting the radial dependence of the gas temperature in each 
of the simulations. The bar in the upper left-hand corner indicates the softening length for 
the low resolution simulation. The horizontal line shows the virial temperature of the cluster 
(Tvir~2.8keV). 
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Fig. 10. — As in Figure 1, presenting the radial dependence of the gas entropy (top) and 
dark matter (bottom) per unit mass in each of the simulations. The bars along the x-axis 
indicate the softening length for the low and high resolution simulations. Note the change 
in scale from logarithmic to linear at IMpc. 

Fig. 11. — A comparison of the gas and dark matter entropy distributions for each of the 
simulations with gas. Each presents the gas (thin line) and dark matter (thick line) entropy 
for Cool+SF, Adiabatic and Lowres simulations. 

Fig. 12. — Each panel plots the projected surface mass density versus radius. The left-hand 
panels plot the total mass in the simulations, while the right-hand panels plot only the dark 
matter. The X, Y and Z's denote the projection axes, which are ahgned with the coordinate 
axes of the simulation. The overdensity is with respect to integrating a uniform density of 
Q = 1 through a line-of-sight of 4Mpc, the size of the box from which these projections were 
made. 

Fig. 13. — The projected X-ray emission: the left-hand panel shows the flux radiated into 
471 steradians in the 2 ^ 10 keV band and the right-hand panel shows the flux in the 
ROSAT/HRI (0.1 ^2.4 keV) band. This band is "windowed" with the normahzed response 
function of the HRI instrument. 
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Fig. 14. — The 2— s>10 keV X-ray surface brightness profiles for the Cool+SF (left) and 
Adiabatic (right) simulations; each panel displays the X, Y and Z cluster projections. The 
thick sohd line in each panel plots the best fit /9-model (eq. 6) over a restricted range from 
50 kpc to 750 kpc, similar to that employed with real X-ray data, while the dot-dashed 
lines plots the best fit /3-model over the full range of the data. To avoid the effects of the 
X-ray peaks in the cluster center and at the positions of galaxies on the mean profile of the 
Coo/-/-5'F simulation, we use only the circled points in the fits. The best fit values for (/3, rc) 
for the Cool+SF simulation are (0.97, 218kpc) and (0.63, 72kpc) for the full and restricted 
range respectively. The corresponding values for the Adiabatic simulation are (1.07, 211kpc) 
and (0.77, 106kpc). 

Fig. 15. — Like Figure 13 but plotting the projected X-ray emission weighted temperature 
profiles. 

Fig. 16. — The projected Sunyaev-Zel'dovich decrement for the various cluster simulations. 

Fig. 17. — The mass accretion rate of gas through spherical shells into the cluster normalized 
by the accretion rate required to accumulate the mass of gas within the cluster virial mass 
at the present time (2.05 x 10^^ M©) in a Hubble time. Note the change of scale from 
logarithmic to hnear at IMpc. 

Fig. 18. — The spherically averaged radial infall velocity of the gas divided by the local 
dynamical velocity, which we define as r/tdyn- Note the change of scale from logarithmic to 
linear at IMpc. 



-59- 



Fig. 19. — The spherically averaged radial infall velocity of the gas in units of the local 
sound speed. Beyond the virial radius a weak shock is present. Note the change of scale 
from logarithmic to linear at IMpc. 



Fig. 20. — The total entropy in the cluster as a function of gas mass between a radius r 
and the virial radius. Both axes have been normalized by the mass of gas within the virial 
radius. 
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Table 1. A summary of previous numerical studies of galaxy clusters. 



Ref. 


Spatial {h ^kpc) 


Dark(^-iM0) 


Virial (/i-iMg) 


Gas(/i-i Mq) 


Cool / Star 


MciuBt(?i"^M0) 


(a) 


150 


3.0x10^^ 


5.6xl0i2 


3.4xl0i2 


N 


1.2 x 10i5(< Ra) 


(b) 


14 


4.8x10" 


4.5 xlO^ 


4.8xlOio 


Y/N 


1.0xlOi4(< il200) 


(c) 


270h p 


l.lxlO"/i 


3.3xl0i3;i3 


2.2x10"/! 


N 


1.85xl0i5;i(< il200) 


(d) 


40 


3.8xl0i2 


l.lxlO" 


4.2x10" 


Y/N 


5.9xl0i4(< R200) 


(e) 


628 


1.2x10" 


4.1x101" 


1.2xl0i° 


N 


many 










(2.1x1012) 






(f) 


78 


2.0x10^2 


7.9x10" 


2.2x10" 


N 


4.7 X 10i'*(< R400) 


(g) 


50 p 


9.4x10" 


2.1x10" 


9.4xl0l° 


N 


2.5 X 10i*(< R200) 


(h) 


50 


4.8x1012 


2.1x10" 


2.6x10" 


N 


3-20 Xl0i4(< Rsoo) 


(i) 


98 


1.6x1010 


1.6x1012 


1.7x10^ 


N 


many 










(3.0x10") 






(j) 


14 p 


2.3xl0i2 


4.5x10^ 


2.6x10" 


N 


1.0xl0i5(< R200) 


(k) 


66 


8.0x10" 


4.8x10" 


4.9x10" 


Y/N 


> 1.2 X 10i*(< Ra) 



Note. All values are for clusters at z=0. Physical, rather than comoving, resolution scales are denoted by a p; corresponding 
virial mass scales as (1 + z)^. The studies are; (a) Evrard (1990); (b) Katz & White (1993); (c) Schindler & Bohringer (1993); 
Schindler & MuUer (1993); (d) Thomas & Couchman (1992); (e) Cen & Ostriker (1994); (f) Metzler & Evrard (1994); (g) Navarro 
et al. (1997); (h) Bartelmann & Steinmetz (1996); (i) Bryan & Norman (1998); (j) Eke et al. (1998); (k) Suginohara & Ostriker 
(1998). 



Table 2. A summary of the various resolutions of the simulations presented in this paper. 



Name 


Spatial {h ^kpc) 


Dark Matter (ft-i M0) 


Virial (h-^Mo) 


Gas (h-iM®) 


Cool/Star 


ndark/ngas/nstar 




Dark 


14 


3.7x10" 


4.6x10'' 


NO 


NO 


18134/0/0 


2400 


Cool+SF 


14 


3.5x10" 


4.6x109 


1.9x10^" 


YES 


18221/11582/6987 


2400 


Adiabatic 


14 


3.5x10" 


4.6 xlO^ 


1.9xl0i° 


NO 


18290/17017/0 


2400 


Lowres 


200 


3.5x10" 


1.3xl0i3 


1.9xl0i° 


NO 


17884/17708/0 


2400 


Hires 


1.4 


1.4x10^° 


4.6x10^ 


NO 


NO 


494486/0/0 


10000 



Cluster Vital Statistics; M^ius = 2.1 x lO^^h"! M©, Rvir = Ih-^Mpc & Vcirc(R.vir) = 1000 kms"! 



Note. — Quantities are defined as in Table 1. Although we quote units to be consistent with Table 1, remember that h = 0.5 and that 
the "spatial resolution" is twice the gravitational softening length. The final column refers to the number of particles, and their type, within the 
virial radius. 
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Table 3. Cluster X-ray properties 



Name Lx(10^3gj.g g-i) (^j^^s^ (^gy) 

Cool+SF 3.14 4.41 

Adiabatic 2.58 3.43 

Lowres 0.88 2.16 



Note. — A summary of the X-ray properties, 
in the 2^10 keV energy band, for the simulated 
clusters presented in this paper. The tempera- 
ture is the mean emission weighted temperature 
averaged over a circle of radius of the virial ra- 
dius (change this to twice or half the virial radius 
does not significantly alter either the X-ray flux 
or emission weighted temperature for the simu- 
lated cluster). The cluster virial temperature is 
2.8 keV. 
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